Maximum-likelihood reconstruction of CP maps 
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We present a method for the determination of the com- 
pletely positive (CP) map describing a physical device based 
on random preparation of the input states, random measure- 
ments at the output, and maximum-likelihood principle. In 
the numerical implementation the constraint of completely 
positivity can be imposed by exploiting the isomorphism be- 
tween linear transformations from Hilbert spaces TL to K, and 
linear operators in K^Ti. The effectiveness of the method is 
shown on the basis of some examples of reconstruction of CP 
maps related to quantum communication channels for qubits. 

PACS Numbers: 03.65.-w, 03.67.-a, 03.67.Hk 



The problem of characterizing a physical device in the 
quantum domain and of reconstructing input-output re- 
lations has been recently addressed in a number of papers 
|H,E|. These issues are obviously interesting for techno- 
logical purposes. For example, the practical determina- 
tion of the transformations acting on quantum states is 
of great relevance in the new fields of quantum informa- 
tion, computation and cryptography H. In this realm 
physical objects as optical fibers, parametric amplifiers, 
directional couplers, C-Not gates, quantum cloning ma- 
chines, quantum communication channels, etc., should 
be characterized with very high precision. From the the- 
oretical point of view a physical device is described in 
terms of a completely positive (CP) map. Hence, an ex- 
perimental method of reconstructing a CP map would 
lead to a complete characterization. On the other hand, 
an effective and reliable technique for the determination 
of a CP map could allow to check experimentally the 
correctness of the theoretical assumptions made in the 
description of the physical device. Finally, recall that 
the structure of CP maps naturally emerges in the the- 
ory of open systems Q . It follows that an experimental 
technique to estimate a CP map also allows to investi- 
gate the interaction between different systems, typically 
a system-reservoir interaction. 

Here in this letter we consider the general problem 
of reconstructing the CP map related to a physical de- 
vice, without any assumptions on its mathematical form. 
We propose a method which resembles ordinary tomogra- 
phy for the use of random quantum measurements at the 
output, but also employs random input states, in order 
to have the richest statistics. The maximum-likelihood 
principle is then applied, using a suitable parameteriza- 
tion of CP maps which is allowed by the isomorphism 
between linear transformations from Hilbert spaces TL to 
/C and linear operators in JC ® TL |g|,|6| . The maximum- 
likelihood method has been used in the context of phase 



measurement [|7|, and to estimate the density matrix |^ 
and some parameters of interest in quantum optics H] . 

The maximum-likelihood principle says that the best 
estimation of some unknown parameters is given by the 
values that are most likely to generate the data one ex- 
perimenter observes. This principle is quantified by a 
maximum search of a functional of the unknown param- 
eters that corresponds to (the logarithm of) the theo- 
retical probability of getting the data one has collected. 
In the following we derive the likelihood functional C{£) 
that links the experimental outcomes with the unknown 
CP map £ that characterizes the device. We consider 
a sequence of K independent measurements on the out- 
put, each described by a POVM Fi{xi), where xi denotes 
the outcome at the /th measurement, and 1 — 1,2, ..., K. 
We denote by pi the state at the input at the /th run. 
The probability of getting the string of outcomes x ~ 
{xi,X2, ■■■,xk} is then given by 



p{x)=Iit{^A^{pi)Fi{xi)] 



(1) 



The maximum-likelihood principle states that the best 
estimate of the map £ maximizes the expression in Eq. 
(0) over the set of completely positive maps. More con- 
veniently, one can search the maximum of the logarithm 
of Eq. (pi), namely 



K 



C{£) =Y,\ogT^[£{pi)Fi{ 



xi] 



(2) 
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The likelihood functional C{£) is concave, and in the 
present case it is defined on a convex set — the set of CP 
maps. It follows that the maximum is achieved by a sin- 
gle CP map, or by a closed convex subset of CP maps. 
In the last case one can infer that the data sample is 
not sufficiently large, or the set of measurements is not a 
quorum. 

The maximization problem is constrained by com- 
pletely positivity and trace-preserving properties of the 
map £. A trace-preserving CP map is a linear map from 
operators in Hilbert space TL [dim(H) = N] to operators 
in /C [dim(/C) = M] which can be written equivalently as 
follows 



£{p)^Y.'^^P'^ 



Tr^KlK®/)^] 



N^ 



= X! ^" ^" p ^» ' 



(3) 
(4) 

(5) 



n=l 



where 

k 

S>0 and Tr^[S] = 1^ , 



N^ 



Tr[l7/(7,]-% and ^p„ t/tf/„ = !„ 



(6) 
(7) 
(8) 



and T denotes the transposition. Eq. (ra) is the weU 
known Kraus decomposition. Eq. (Q) exploits the iso- 
morphism between hnear maps from Ti. to IC and hnear 
operators on the tensor-product space K. <^H |5|,P,p|,|lO| . 
The operator S can be written in terms of the map £ as 
follows 



N^ 



^ = ^5(i/,)(^y; 



(9) 



where * denotes the complex conjugation, and {Vi} is 
any orthonormal basis for the space of linear operators 
on Ti, namely 



and for any operator O 



N^ 



= ^Tr[y;0]l-, . 
Notice also that for linearity ||9|,|ll[] 



(10) 



(11) 



(12) 



where |^) represents the (unnormalized) maximally en- 
tangled state |\1/) = X]n=i l*^)!")- 
Eq. (0) can also be written as 

£{p)=Tin[{^K(^p)S^], (13) 

where T denotes the partial transposition in Ti, and then 



5^ = ^ £(!/.) ®V;t 



(14) 



Finally, Eq. (||) can be shown as follows. Chosen an 
orthonormal basis {Vi} in the sense of Eqs. (10) and 



N^ 



(y), Eq. (g) rewrites £{p) = J2i,j=i 'l^jV^pVJ , with q,j = 
J2k Tr[^feT^/]Tr[^|.l^j]. The matrix Q with element q^j is 
positive, and it can be diagonalized Q — WDW'^ with 
W unitary. Then one has 

£ (P) = E {WDW% V, PV^^Y. P" Ur^ P ^n , (15) 

i,j — l n— 1 

where p„ — {D)nn > and {Un} is the new orthonormal 
basis Un = Y.iLi{W)ikVi. 



For Ti, = JC the matrices Ak and f7„ are squared. When 
referring to quantum communication channels, the chan- 
nel is called bistochastic if also J2k ^kAl — tn- More- 
over, operators {Un} in Eq. dq) could be unitary, and in 
such case the (bistochastic) channel is said to be given 
by external random fields. For a qubit system (Ti — C^) 
the set of bistochastic and external-random-field channels 
coincide III4 . 

The isomorphism between linear maps and operators 
established in Eqs. (|)-(|) or ( |l3[)- (|il) has been use- 
ful for the study of positive maps | lEOL and to address 
the problem of separability of CP maps [Q. For our 
task — the maximization of C{£) — it is also crucial. The 
condition S > allows one to write 



s = C^C 



(16) 



where C is an upper triangular matrix. Moreover, the 
diagonal elements of C can be chosen as positive. Such 
decomposition — referred to as Cholesky decomposition — 
is commonly used in linear programming |l3]. Similarly, 
one has for the matrices pf and the POVM's Fi{xi) 



Pi 



R]Ri 



Fi{xi) ^ A\[xi)Ai(xi) 



(17) 



From Eqs. 
tional in Eq. 



( [iq ) and (p7|), the log- likelihood func- 
I) rewrites 



K 



C{£) = C{C) ^Y.^ogTT[C^C{R\Ri ® A\{xi)Ai{xi))] 



K 



1 = 1 
NM 



5] log E \in\C{R\®A\{xi))\m)) 



/^l n,m— 1 



(18) 



where {|?^))} denotes an orthonormal basis for Ti ® JC. 
The expression obtained in Eq. (HSh for the likelihood 
functional automatically satisfies the constraint of com- 
pletely positivity. Furthermore, the terms appearing as 
argument of the logarithm are explicitly positive, thus 
assuring the stability of numerical methods to evaluate 

The trace-preserving condition is given in terms of the 
matrix S as Tr;c:['5'] = !-«. This can be taken into ac- 
count by using the method of undetermined Lagrange 
multipliers, then maximizing 



N 



£'(£) = £(£)- E My (*|Tr^[^]b-) 
= C{£) - Tr[(lK ®p)S], 



(19) 



^N 



where p is the undetermined matrix p = ^^ ■ -^ Pij\i){j\- 
The multipliers pij cannot be easily inferred, except the 
condition Tr[/x] = K . Writing S in terms of its eigen- 



vectors as S* = ^is1\si 



, the maximum likelihood 



condition dC'{£)/dsi = implies 



K 



y;pTT[{pf(g,Fi{xi))s,\s^ 



1=1 



Tiiipf (g> Fiixi))S] 



Tr[(lK (g) p)s, \s. 



(20) 



Multiplying by Si and summing over i gives Tr[/i] = K . 
However, notice that the constraint Tr[S'] = N which fol- 
lows from TrK[S'] — tn isolates a closed convex subset of 
the set of positive matrices. Hence, the maximum of the 
concave likelihood functional still remains unique under 
this looser constraint, and one can check a posteriori that 
the condition TrK;[5'] = 1-h is fulfilled. The functional we 
maximize is then 



C{C)=CiC)-^Tr[C^C] 



(21) 



where C{C) is given in Eq. ( |18| ) , and the value of the mul- 
tiplier has been obtained through a derivation similar to 
Eq. (pQ). The number of unknown real parameters is 
given by (NM)^. The problem of maximization of func- 
tional as Eq. (UW enters the realm of programming and 
numerical algebra optimisation, where various techniques 
are known ]14| . 

In the following we show the effectiveness of our 
method on the basis of some examples of quantum com- 
munication channels for qubits, i.e. Ti ~ K, = C^ . In this 
case {NMY — 16, and for such a relatively small num- 
ber of parameters one can efficiently apply the method 
of downhill simplex |l4,n3] to find the maximum of the 
log-likelihood functional. This method has been reliable 
in the reconstruction of the density matrix of radiation 
field and spin systems P|. The results in the following 
simulations have been obtained using random pure states 
at the input of the channel, along with a projective mea- 
surement in a random direction at the output. 

The first example is the Pauli channel for qubits 



^p(p) = zL'p^'^'^p^'- ' 



(22) 



where X]i=i-P« ^ ^, o-q = Ij and Ui (i — l^jS) denote 
the customary Pauli matrices. From Eq. (|2|) the corre- 
sponding positive matrix Sp writes 



Dn 



Po+Ps P0~P3 

P1+P2 P1-P2 

P1-P2 P1+P2 

, Po - P3 P0+P3 



(23) 



on the lexicographically ordered basis |00), |01), |10), 
1 11), where |0) and |1) corresponds to the eigenstates 
of az with eigenvalues 1 and —1, respectively. 

In table | we reported the reconstructed matrix ele- 
ments of Sp as obtained by a Monte Carlo simulation with 
K = 30000 runs, for theoretical values po = 0.3, pi = 0.2, 
P2 = 0.4, and ps — 0.1. The trace-preserving property 
corresponds to the conditions 5*^(1,1) -I- 5*^(3,3) = 1, 
Sp{2, 2) + SpiA,4) = 1, and 5*^(1, 2) -t- 5*^(3, 4) = 0, which 
are clearly satisfied. The estimated values compare very 
well with the theoretical ones. 

For pi = p2 — P3 in Eq. ( p^ ) , one obtains the depolar- 
izing channel 



fd(p) = Ap+^— 1 



(24) 



with A = 1 — 4pi. In Fig. |l] (circles) we reported the sta- 
tistical error {5X)ml in the evaluation of the parameter A 
versus the size K of data sample (with theoretical value 
A = 0.8). The value of Xml has been inferred by the 
complete reconstruction of the matrix Sp. However, no- 
tice that one can also implement the maximum-likelihood 
method upon assuming the form of the CP map as in 
Eqs. (P^ ) or (gj). In such case the space of parameters 
is reduced to 4 and 1, respectively. Fig. Il| also shows the 
results obtained by a 4-parameters estimation (triangles), 
thus by assuming an external-random-field channel. In 
both cases one has an asymptotic inverse square-root de- 
pendence of the statistical error on the size of the data 
sample ((5A)j\/l oc K^^^^, according to the central limit 
theorem. 
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FIG. 1. Maximum-likelihood estimation of the CP map re- 
lated to the depolarizing channel. The picture shows the value 
of the statistical error {S\)ml in the estimation of the pa- 
rameter A (theoretical value A = 0.8) versus the size K of the 
data sample. Circles referred to a ML reconstruction without 
assumptions on the form of the CP map; triangles are the 
results when assuming the external-random-field form. The 
asymptotic dependence of the statistical error versus K is in- 
verse square-root {5X)ml oc K~^'^, as it is demanded by the 
central limit theorem. 

In the last example we consider a non-bistochastic 
channel, namely the amplitude damping channel 



£aip) ^ Ml p Ml + M2 p M2 



with 
Ml 



1 

Vp 



and M, = 



VT^ 





(25) 



(26) 



The corresponding positive matrix Sa write 



Sa 



1 0^' 

l~p 


Vp p 



(27) 



In Fig. 1^ we have plotted the estimated value pml of pa- 
rameter p its theoretical value, as inferred by the recon- 
struction of the matrix Sa through K = 10000 random 
measurements. 



FIG. 2. Maximum-likelihood reconstruction of the CP map 
describing an amplitude damping channel. The value of the 
parameter p is inferred by reconstructing the positive matrix 
Sa in Eq. (p7|). Random pure states at the input and projec- 
tion along random directions at the output have been used, 
with K = 10000 measurements. 

In conclusion, we have proposed a method for recon- 
structing the completely positive map related to a phys- 
ical device, based on the maximum likelihood principle. 
The method is very general, does not require a priori 
knowledge of the mathematical structure of the CP map, 
and can be adopted in many fields as quantum optics, 
spin systems, optical lattices, atoms, etc. We have shown 
some examples of reconstruction of CP maps related to 
quantum communication channels, applying the down- 
hill simplex method for the search of the maximum of 
the likelihood functional. 
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TABLE I. Maximum-likelihood estimation of real and 
imaginary parts of the matrix elements of Sp related to the 
Pauli channel in Eq. (p2), for theoretical values po = 0.3, 
pi = 0.2, p2 = 0.4, and pa — 0.1. Random pure states and 
projective measurements along random directions have been 
used, with K = 30000 runs. Compare with the theoretical 
values as obtained from Eq. (Eq). The statistical error in the 
estimation of the matrix elements is around 0.01. For typical 
values and behavior of the statistical errors, see Fig. nl 



